Cluster the data
fviz_nbclust(
table_cluster_scaled,
pam,
k.max = 10,
method = "wss",
diss = get_dist(table_cluster, method = "spearman")
#diss = get_dist(table_cluster, method = "euclidean")
)

n_clust = 6
#k_n <- kmeans(table_cluster_scaled, centers = n_clust, nstart = 100, iter.max = 15000)
k_n = cluster::pam(table_cluster_scaled, k = n_clust)
fig_cluster <- fviz_cluster(k_n, data = table_cluster_scaled, ellipse = T, palette = "Set2",
geom = "point") +
theme_classic()
fig_cluster

df_cluster_info <- tibble(SAMPLE_ID = row.names(table_cluster_scaled),
cluster_id = k_n$cluster)
bangladesh_gw <- bangladesh_gw %>%
left_join(df_cluster_info)
## Joining with `by = join_by(SAMPLE_ID)`
bangladesh_gw <- bangladesh_gw %>%
mutate(cluster_id = factor(cluster_id))
bangladesh_gw %>%
filter(!is.na(cluster_id)) %>%
ggplot(aes(x = cluster_id, y = As_ugL, fill = cluster_id)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set2") +
scale_y_log10() +
theme_classic() +
theme(legend.position = "none")

bangladesh_gw %>%
filter(!is.na(cluster_id)) %>%
ggplot(aes(x = cluster_id, y = Fe_mgL, fill = cluster_id)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set2") +
scale_y_log10() +
theme_classic() +
theme(legend.position = "none")

bangladesh_gw %>%
filter(!is.na(cluster_id)) %>%
ggplot(aes(x = cluster_id, y = SO4_mgL, fill = cluster_id)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set2") +
scale_y_log10() +
theme_classic() +
theme(legend.position = "none")

bangladesh_gw %>%
filter(!is.na(cluster_id)) %>%
ggplot(aes(x = cluster_id, y = Si_mgL, fill = cluster_id)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set2") +
scale_y_log10() +
theme_classic() +
theme(legend.position = "none")

bangladesh_gw %>%
filter(!is.na(cluster_id)) %>%
ggplot(aes(x = cluster_id, y = WELL_DEPTH_m, fill = cluster_id)) +
geom_boxplot() +
scale_fill_brewer(type = "qual", palette = "Set2") +
scale_y_log10() +
theme_classic() +
theme(legend.position = "none")
## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

summary(bangladesh_gw$cluster_id)
## 1 2 3 4 5 6 NA's
## 356 599 731 639 811 393 5
sf_bangladesh_gw <- bangladesh_gw %>%
st_as_sf(coords = c("LONG_DEG", "LAT_DEG"))
tmap_mode("view")
## tmap mode set to interactive viewing
sf_bangladesh_gw %>%
filter(!is.na(cluster_id)) %>%
tm_shape() +
tm_dots(col = "cluster_id", palette = "Set2",
popup.vars = c("cluster_id", "WELL_DEPTH_m", "As_ugL")
) +
tm_mouse_coordinates() +
tm_basemap(server = c("Esri.WorldImagery", "OpenStreetMap", "Esri.WorldShadedRelief")) +
tm_scale_bar()
## Warning: Currect projection of shape . unknown. Long-lat (WGS84) is assumed.
sf_bangladesh_gw %>%
filter(!is.na(cluster_id)) %>%
tm_shape() +
tm_dots(col = "As_ugL", palette = "viridis", breaks = c(0,5,10,50,100,200),
popup.vars = c("cluster_id", "WELL_DEPTH_m", "As_ugL")
) +
tm_mouse_coordinates() +
tm_basemap(server = c("Esri.WorldImagery", "OpenStreetMap", "Esri.WorldShadedRelief")) +
tm_scale_bar() +
tm_facets(by = "cluster_id", sync = T)
## Warning: Currect projection of shape . unknown. Long-lat (WGS84) is assumed.
## Warning: Values have found that are higher than the highest break